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Abstract. We study the one-dimensional diffusion process which takes place between 
two reflecting boundaries and which is acted upon by a time-dependent and spatially- 
constant force. The assumed force possesses both the harmonically oscillating and the 
constant component. Using techniques presented in the paper 1^ we derive the set 
of integral equations whose comprise the Green function for our diffusion process. In 
the later part we focus on the time-asymptotic stationary regime of considered non- 
equilibrium process. Some of its time-asymptotic characteristics, e.g. time-averaged 
and time-asymptotic probability density of the particle coordinate, exhibit nontrivial 
features. 



1. Introduction 

Consider a rectangular vessel filled with small Brownian particles dissolved in a good 
solvent. With no external force acting on our particles, the concentration of the particles 
would be homogeneous. Now assume that each particle carries the same electric charge. 
What would be the concentration profile if we put the vessel in an external, constant and 
homogeneous electric filed? After some time it approaches the new constant equilibrium 
profile. Its shape would depend on many factors, like temperature of the solvent, 
dimensions of the vessel, mobility of the particles and a strength of the field. And 
what if the intensity of the filed would oscillate? 

We assume that the external homogeneous field acts only along the one axis of 
the vessel, let say along the x-coordinate, that the particles do not interact between 
themselves and that the concentration of the particles is homogeneous in the other {y 
and z) directions. Under these circumstances we can treat the emerging dynamics as 
a one- dimensional diffusion process confined to the spatial interval by two reflecting 
boundaries and acted upon by a space-homogeneous and time-dependent force. 

Since the beginning of the 20*'^century the diffusion problems concerning the 
movement of the Brownian particle in a fixed potential profile have been investigated 
many times in the literature and nowadays it is a well understood problem |21 El E] • In 
X corresponding author 
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a past few decades the diffusion dynamics was reexamined in the systems, in which the 
potential incorporates a time-dependent component. Notably interesting phenomena 
such as a stochastic resonance [HI E] or effect of Brownian motor [Zj arise if the potential 
changes periodically in time. It is known experimentally that nonlinear effects of 
alternating electric field are very important in governing the function of membrane 
proteins [SI E] ■ 

The time-dependency of the potential represents the fundamental mathematical 
difficulty in the theoretical treatment of the diffusion problem. Nowadays there are only 
few exactly solvable diffusion problems with the time- dependent potential [Tm[THIT^ I13j. 
Moreover, the presence of a (reflecting/absorbing) boundary rises further the complexity 
of the computations. 

The next section presents the general treatment of the problem based on the 
computational techniques presented in p. In section 3 we assume that the force, aside 
from its constant part, incorporates the harmonically oscillating component. Section 4 
presents the procedure, which allows the improved discussion of the time-asymptotic 
dynamics. Section 5 is focused on the time and space resolved probability density 
and its period-averaged counterpart. Section 6 contains the concluding remarks and 
summarizes the results. 



2. Green function for a diffusion between two boundaries 



Consider a free diffusion of a Brownian particle under the influence of a space- 
homogeneous and time- dependent force F(t). Such setting corresponds to the time- 
dependent potential of the form V{x,t) = —xF{t). In the overdamped regime, the 
time-evolution of the probability density for the position of the diffusing particle is 
determined by the Smoluchowski equation jl] 



d , , d 



d 

-D—p{x, t) + v{t)p{x, t) 



(1) 



Here, v{t) is the time-dependent velocity of the drift motion, which is induced by the 
external force, v{t) = F{t)/T. Parameter F equals the particle mass times the viscous 
friction coefficient. Its reciprocal value plays the role of the particle mobility. The 
thermal-noise strength parameter D is the linear function of the absolute temperature, 
D = k^T/r. The expression in the square brackets on the right-hand side of (0) 
represents the probability current. 

The solution of Smoluchowski equation ((H) with the assumed time- dependent 
potential is accessible using the various algebraic methods O US] reads 

2 



G{x,t;x',t') 



1 



^ADTx{t-t'] 



exp 



x-x' - /; dt"v{t") 



4D{t - 1' 



(2) 



This Green function yields the solution of (Q) for the initial condition at time t = t' in 
the form limt^f p{x,t) = 6{x — x'). Qualitatively, equation (j2I) represents the gradually 



Diffusion on the interval 



3 



spreading Gauss curve, whose center is located at the coordinate x = x' + j^, dt"v{t"). 

Now assume, that there are two reflecting boundaries located at Xi and X2, with 
Xi < X2- Let the particle is initially placed at the arbitrary coordinate between the 
boundaries, i.e. x' G {x\,X2)- Without loss of generality we set the time origin at t = 0. 
The resulting diffusion dynamics is restricted by the boundaries to the space-interval of 
length I = X2 — X\. For the Green function of the present problem we assume the form 

?7(x,t;x',0) = G(x,t;x',0) - D \ dt' —G{x,t; xi,t')Biit' , x') + 



D / dt'—G{x,t;xi,t')Bi{t',x') 
Jo 

D dt'^G{x,t;x2,t')B2{t',x') , (3) 



+ 

for all x,x' G {xi,X2) and t > 0. Mathematically, our diffusion problem reduces to the 
determination of two unknown functions Bi{t,x') and B2{t,x'). 

The next steps follow the procedure that was presented in great mathematical 
details in PJ. The reader interested in the mathematical subtleties is referred to the 
mentioned paper. The requirement of vanishing probability current through every 
reflecting boundary, together with our ansatz 0, leads to the set of two Volterra integral 
equations of the flrst kind for the unknown functions Bi{t,x') and B2{t,x'), 



t 

Dl dt' 

'0 



G{xi,t; Xi,t') 
-G{x2,t;xi,t') 



-G{x2,t;xi,t') 
G{x2,t;x2,t') 





' Bi{t',x') ' 






B2{t',x') 





P^dxG{x,t;x',0) 
tdxGix,t;x',0) 



(4) 



It can be shown, that both solutions have a speciflc physical meaning. Each one 
represents the actual probability density at the respective boundary, i.e. Bi(t,x') = 
U{xi,t; x', 0) and E2(t, x') = U{x2, t; x', 0). 

The solutions of the system yields the complete description of the emerging 
diffusion dynamics. From them one can from in derive all one-time characteristics of 
our model. This general result holds for an arbitrary time dependency of the potential. 



3. Harmonically oscillating force 

In this preparatory section we concertize the general results of the preceding paragraphs 
to the speciflc time-protocol of the potential force F(t). We assume that it oscillates 
harmonically with a time-independent amplitude Fi > and frequency u around a flxed 
value Fo > 0. The cases with Fq < and/or Fi < are easily accessible from the model 
symmetries. Hence the force induced drift velocity is prescribed by the formula 

v(t) = f + f 1 sm{u!t) , (5) 

where Vq = Fq/T and Vi = Fi/T respectively. To make the following calculations more 
transparent, we introduce the substitution r 

Vo l_ 

2D ^/4^ 



\vo\H/{4:D). The formula turns to 



G{x,t;x',t') 
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(6) 
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with 9 = 4:Du/vq, k = vqVi/ {2Duj), A = VqI/{2D), ^ = x/l and I = X2 — Xi. Here we 
assume uj > and D > 0. 

For a moment, we focus on the case of time-independent force, i.e. we put f 1 = 
in (0). Integrals in and Q now depend on time only through the difference t — t' . 
Therefore the complete solution can be readily obtained using the Laplace-transform 
method or using the expansion into the eigenfunctions of the evolution equation |16j . 
The probability density [/(x, t; x', 0) then relaxes to the equilibrium distribution 

rr{eq)( N ^0 ^Xp [- ^ (^2 - x)] r ^/ \ (r,^ 

^ = D [1 - exp(-2A)] ' " ^ ^"^'"^^ • 

Through this formula, the asymptotic values of the probability density at the boundaries 
B^'^'' and b'^''\ and hence the asymptotic solutions of the system Q are readily 
accessible. Specifically 



A 



5i-) = ^^^ and 5^^ = ;?^^^. (8) 

These formulas can be used to obtain the solution in the term of an adiabatic 
approximation [El HHj. Simply changing vq to v{t) = vq + sin(u;t) and A to X(t) = 
v{t)l/{2D) in (jH)) the adiabatic solutions B[''^\t) and B^2'^\t) are obtained. However, 
such approximation is suitable only for high temperatures, low frequencies and weak 
driving force Fi, i.e. if the system can nearly instantaneously relaxes to the equilibrium 
state corresponding to the actual potential profile. 

Now we return to the original setting Fi 7^ and 7^ 0. In this case, the equilibrium 
distribution U^^'^\x) is never reached. The system asymptotically approaches the non- 
equilibrium steady state. The system integral equations (jH) now represents the principal 
analytical difficulty. The integrals have no more the convolution structure and the 
kernels display the weak singularity at the upper integration limit. Regardless we 
have developed analytical methods, which allow us to make some statements about 
the solution. These results will be presented in the later sections. 

For the numerical solution of (jlj) one can use the product integration method 
[19, 2011^, which treats properly the weakly singular kernels. Its outcome incorporates 
both the transitory effect and the stationary regime. Since we are interested only in 
the asymptotic steady-state and we want to have the better analytical insight into the 
long-time behavior, we use another approach. 



4. Long— time limit of the densities at the boundaries 

In this section we present the analytical method for solving the system of integral 
equations (jH). Though we cannot yet write the exact expression for its solutions Bi(t, x'), 
i = 1,2, our method allows us to gather more information of the solutions, mainly about 
its time-asymptotic behavior. 

The system of integral equations represents the fundamental difficulty in the 
analytical treatment of the emerging diffusion dynamics. The non-convolution structure 
of the equations rules out a direct application of the Laplace transformation. Another 



Diffusion on the interval 



5 



/dwexp | — (u^ — 2in)(r — r') — mK[cos(6'r) — cos(^^r')] — iu—^{x — x')\ 



difficulty represents the second power of the nominator in the exponential which appears 
in (j2I). We got rid of it only at the cost of introducing the Fourier transformation into 
our calculation. Specifically, we use the integral expression for the function (jH]) 

G{x,t;x',t') = (9) 
AirD 

As a next step we expand the exponentials of the cosine into the sums of modified 
Bessel functions Ifc(a), 

oo oo 
g-mKcos(er) _ Ifc(-iMK)e''^^" , ^C^-') = ^ Ifc(iM/t)e''=^"' . (10) 

k=—oo k=—oo 

Now we turn our attention to the right-hand sides of the integral equations. The 
right-hand side of the ffist equation in (jH), let designate it ri(t,x'), represents the total 
probability of finding the particle left from xi, i.e. left from the left boundary, in the 
diffusion problem with the same potential but without any boundary. Analogously, 
the right-hand side of the second equation, r2{x',t), represents, the total probability 
of finding the particle to the right from the right boundary placed at X2- In out case 
with no boundary and f o > the particle gradually tends to the plus infinity. Hence 
asymptotically the particle will be surely found the right from the considered interval. In 
other words, the integral on the right-hand side of the ffist equation in Q asymptotically 
vanishes. At the same time the right-hand side of the second equation approaches the 
unity. According to this observation, we expect the right-hand sides to have a form 

ri(t, x') = fi(t, x') , r2(x', t) = l + r2{t, x') . (11) 

where rj(t) are the transitory parts, that have the important property lim^^oo 'f^i{x', t) = 
0, i = 1, 2, for any x'. 

From the assumed form of the right-hand sides, we induce the ansatz for the 
solutions Bi(t,x'), i = 1,2, i.e. we assume 

B.,it,x') = ^ [/5f)(t) + A(t,x')] , 1 = 1,2. (12) 

Again, we suppose that the transitory parts (3i{t,x') vanish identically for t — » oo. 
Moreover, from the assumption, that the asymptotic solutions are the periodic function 
of time with the fundamental frequency uj, we expect the asymptotic parts to have a 
form 

oo 

= P^,kexp{-ikut) , 1 = 1,2. (13) 

A;=— oo 

Here (3i^k are the unknown complex numbers and we refer to them as to the complex 
amplitudes. 

After all above described preparative steps, we can treat the original system (0]) 
with the Laplace transformation. Using the fact, that the exponentials from (|Tn|) and 
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m,n=— 00 ^ ^ 

IP^^ g I^(-it..)I (i^.)e-^ (,)^^^^^^ ^ ^ ^^^^ 

m,n=— 00 I \ /JL \ /J 

1 /-"^ ^ I^(-im)I„(iM/t) X , 1 -/ 



TT 



m,n=—oo 



Here we introduced the abbreviations = ±i[vT~-i-^m =t 1] and = z — imO. 

In the above equations we now proceed with the w-integration using the residue 
theorem. For the equation (|14p the integration path has to be closed in the upper half 
of the complex plain. The only singularities there lies in u+{zm+n)- Contrary to this, 
the integration paths in ()15|) must be closed in the lower complex half-plane and the 
relevant singularities are located at u_{zm+n)- 

Now we carry out the time-asymptotic limit in the 2;-domain. We multiply both 
sides of the equations by the Laplace variable z and perform the limit z —>■ 0^. 
According to our assumptions, the functions ri{z,x') and l3i{z,x'), i = 1,2, have no 
pole at the origin and they vanish in the considered limit. The only other pole at 
2; = is incorporated in the function f3^'^\zm+n) and according to (fTSj) we can write 
lim^^o+ z(3^i''\zm+n) = A,m+n for i = 1,2. 

The resulting equations read 

m,n=—oo * m,n=— oo * 
- e ™ . /?l,m+n+ > . . „ /:/2,m+n = 1 • (17) 

^ — ' Vl — ifn" — V 1 — I'^o' 

m,n=—oo * m,n=—oo * 

Here = a/1 — im^±l. The above rather cumbersome expressions can be transformed 
into the compact matrix form 

LiMoMi - LiMiRi |/?2) = |0) , (18) 

-L2M2M2 + L2M0M2 1/92) = 1 |0) . (19) 

The matrix elements of Lj and are build up from the Bessel functions, 

(m| Li \n) = l\rn-n\ (f^an) , ("^1 L2 |n) = l\rn-n\ (-Z^"^) , (20) 

(m| Ml \n) = I|„_„|(-/ta+) , (m| M2 \n) = l\m-n\ifio:;^) ■ (21) 
and the diagonal matrices Mj are defined as 

(m|Mo|n) = ^=^^ , (rr,|M,|^) = (m|M2|n) = ^^^^== . (22) 

V 1 — im^ V 1 — im9 V 1 — im^ 
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The detailed analysis of elements of matrices Lj and Mj reveals the possibility to 
rearrange (j^^ into the simpler form 

Mi|/3i) -DiMil/32) = 0|0) , (23) 

-D2M2 \f3i) + M2 1/^2) = 1 |0) . (24) 

with 

(m| Di \n) = 5mn e"^"™ , (m| D2 \n) = 6^n e"^"^ • (25) 

At present time we cannot write the explicit solutions {^2) of the system (j^ - P^ . 
Is shows up, that the complex amplitudes Pi^k and Pi-k are complex conjugated and 
is always real and positive. The solution can be written in the form 

Pi''\t) = A,o + 2 V P^,k cosikcot - 0,,fc) , = arctan f ^^M) ' ^ = 1.2. (26) 

Hence the asymptotic solution of the original system of integral equations (j^ are 
Bl^\t) = %Pt\t),t = l,2. 

The formula directly imply, that the period-averaged time-asymptotic density 
at the left or the right boundary is b["'^^ = ^/3i,o or 32""^^ = ^^2,0 respectively. It 
can be proved that there is a tight connection between these values. After multiplying 
the both equations of the system (j^ - P^ by the vector (0| from the left and after 
evoking the properties of matrices involved, one equation reduces to the simple statement 
^2,0 — Pi,o = 1- Differently speaking, the difference between time-asymptotic and time- 
averaged values of the density at the right and the left boundary always equals to the 
ratio vq/D. 

The fundamental difficulty in the calculation of the complex amplitudes is the 
complicated inversion of the matrices Mi and IR2- We have solved the system (j23|l - 

numerically and inserted the outcome into (f^H|) . The results are in the complete 
agreement with the later-time outcome of the product integration method. 



5. Time-asymptotic probability density 

Assume the vectors of complex amplitudes and {^2) are known. Now we turn out 
attention to the asymptotic behavior of the whole probability density, i.e. we focus on 
the function U{x, t; x', 0) ~ U^'^^x, t) for t — > 00. In a tight analogy with the ansatz for 
the densities at the boundaries (fT^ we expect U^'^^ {x, t) in the form of its Fourier series 

+00 

U^''\x,t) = ^ Yl M^) (^wHkujt) . (27) 
fc=— 00 

Of course, it must hold that U^^'^xut) = B^^\t) and U^''\x2,t) = Bi''\t). To 
determinate the x-dependent Fourier coefficients in (jTTj) we use the similar "tricks" 
as for the calculation of complex amplitudes. We start with the equation Q. The 
first term on the right hand side, G{x,t; x' ,0), for a fixed x represents the gradually 
diminishing oscillations and vanishes as t — 00. In the other two terms, we use the 
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X [m] 

Figure 1. The time and space resolved probability density in the asymptotic regime 
plotted over two periods of external field. Parameters used are vq = 0.1 ms~^, 
vi = 4.0ms~\ Z = 4.0m, = 4.0rads~^ and D = l.Om^s"^. 



expression Q for G{x,t; Xi,t') and proceed with the x-derivation. From this point the 
procedure is the same as in the preceding section, only the proper integration contours in 
the complex plain must be chosen during the ^-integration. Skipping the computational 
details, the Fourier coefficients Vk{x) = {k\v) are given by the formula 

It;) = L2N1M2 lA) + L1N2M1 I/32) . (28) 

The x-dependence of coefficients stems from the matrices Nj. Their elements are given 
by the expressions 



(m| Ni \n) 
(m| N2 \n) ■■ 



1 



^/l — imO 



- 1 



2 [VI -im^ 



exp 



exp 



— Xi)a^ 



'2D 



{x2 — x)a^ 



(29) 
(30) 



The matrices Lj, Mj and the combinations are introduced in the preceding section. 
Hence if we know the complex amplitudes, i.e. the vectors and {(^2), we can from the 
expression calculate the coefficients Vk{x). The result of our numerical calculations 
is displayed in the Figured 
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01234 01234 
X [m] X [m] 

Figure 2. The comparison of time-asymptotic probability density averaged over the 
period of the external field U'^°''"\x) (full line) and the time-asymptotic probability 
density in the static case U^'^'^\x) (dashed line). Parameters used for the calculation 
of the curves are in the left panel the same as used in figure ^ vq — O.fms^^, 
vi = 4.0ms~"'^, / — 4.0m, u — 4.0rads~"'^ and D = l.Om^s"^. The only parameter 
changed on the right panel is the driving frequency uj — 2.0rads~^. 



Knowing the time-resolved probability density one could ask, how does look the 
average density profile. If we perform the time-averaging directly in (j2Z|), we readily 
recognize, that the average probability density 1/^°''"^ (x) corresponds to the zeroth Fourier 
coefficient vo{x). 

jjiav) = ^\im " dt' U{x, t'; x', 0) = Vo{x) . (31) 

The Figure El demonstrates the numerically calculated averaged probability density 
profile U^'^^ (x) and compares it with the time-asymptotic probability density U^'^'^^ (x) in 
the static C3iSG clS defined in (jZj). For certain sets of model parameters different shapes 
of U^'^^x) can occur. Two of them are shown at Figure |21 We discuss them in the next 
section. 



6. Concluding remarks 

Let now think again about the vessel in the external oscillating field that we have 
outlined in the beginning of the paper. We have presented rather complete picture of 
the long-time dynamics of the concentration profile. The solutions of the system of 
integral equations Bi{t,x') and Bi{t,x') (or their asymptotic counterparts B['^\t) and 
B^\t)) can be interpreted as a relative portion of the particles touching the respective 
wall of the vessel. 

The space and time resolved concentration profile in our vessel can be possibly 
measured, e.g. by measuring the light absorption or transmission in the direction 
perpendicular to the x-axis. If the measuring device with higher time-resolution than 
the driving frequency u would be used, one will obtain the results analogous to Figure UJ 
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i.e. the function U^^\x, t). On the other hand, if the time- resolution of the device would 
be several or more cj, one can measure only the averaged concentration, i.e. the time- 
averaged density profile f/("^) (x) demonstrated at Figure El 

Our numerical calculation has revealed the interesting features in the behavior of 
the time-averaged asymptotic density U^"''"\x). In the non-static mode, i.e. wi 7^ and 
u; 7^ 0, the average probability density in the vicinity of the boundaries is always higher 
than in the static case. More precisely, U^'''"\xi) > U^'^i^Xi), i = 1,2. Differently 
speaking, the oscillating external field always rises the concentration at walls of the 
vessel. The equal sign holds in the limit I ^ 00. This case is equivalent to the model 
discussed in p. 

The numerical results show two distinguishable shapes of the function (x). The 
realization of the particular shape is a consequence of the fine interplay of the model 
parameters. Both possibilities are demonstrated at Figure |21 The left panel shows the 
surprising situation when two local minims surround a local maximum and 
crosses [/(^^^ four times. This case can only realize when the oscillating component of 
the field reverses the global slope of the potential during some part of its period, i.e. 
Vo < vi and simultaneously frequency u is sufficiently high and the spatial interval 
between the boundaries is sufficiently long. In comparison with the static case, the 
additional oscillating field rises the average concentration in the middle-region and near 
the walls of the vessel and creates two regions with the lowered concentration. 

The right panel of Figure |21 demonstrates the situation when the driving field acts 
as a centrifuge. The function U^"''"'^ (x) crosses U'-'"^\x) twice. The f/^'^'') (x) can be both 
monotonically decreasing or forming the global minimum in the middle-region of the 
interval. The average effect of the oscillating field is the increased concentration near 
the walls of the vessel and diluted region in nearby its middle. Again, the realization of 
a particular case depends on a tuning of the model parameters. 

The more detailed discussion of the asymptotic behavior of the system and of its 
time-averaged characteristics is obstructed by the rather complicated inversion of the 
matrices Mj, i = 1,2. It is essential during the calculation of the complex amplitudes. 

The presented mathematical procedure is not suitable only for the calculating the 
time resolved or time averaged probability density. Calculating the complex amplitudes, 
i.e. the asymptotic solution of the original system of integral equations, and going 
through analogous presented in Sections 4 and 5, it is also possible to check the proper 
normalization of U^"''^{x, t) or to calculate the mean coordinate. Our present paper deals 
with the kinetic characteristic of the model. Another possibility would be to study its 
thermodynamics of this isothermal non-equilibrium process. 
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